#------------------------------------------------------------------------------
# Import Library
#------------------------------------------------------------------------------

rm(list=ls())
library(LalRUtils)
libreq(data.table, fixest, ggplot2, patchwork, tictoc, RPushbullet, 
       foreach, doParallel)
set.seed(42)
theme_set(lal_plot_theme())
notif = function(x) pbPost("note", x)
# Use cluster to run this code since it takes too much time
#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )


#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# VCF data
#------------------------------------------------------------------------------

# %% subset to primary sample
vcf_data = vcf_data[year>=1990]
vcf_data[, t := year - 1990]
ex_ante_med = quantile(vcf_data$cover_1990, 0.5)
above_med = vcf_data[cover_1990 > ex_ante_med]
above_med[, never_treated := max(D) == 0, cellid]

# %%  timing placebo
# extract upper bound for random timing
ub_by_cell = above_med[, .(ub = first_pesa_exposure[1]), cellid]
# %% function to randomly draw treatment timing in pre-period, re-estimate
permuteTiming = \(){
  # draw random year between 1990 and actual treatment date to switch on treatment
  ub_by_cell[, pseudo_treat := Map(\(ub) sample(1990:ub, 1), ub), "cellid"]
  # merge it back and compute treatment effect
  dftemp = merge(above_med, ub_by_cell, by = 'cellid')
  # fake treatment timing, real scheduled status
  dftemp[, Dfake := (year >= pseudo_treat) * sch]
  m2 = feols(forest_index ~ Dfake | cellid[t] + styear, data = dftemp, cluster = "blk")
  m2$coefficients
}


# Set the number of cores to be used for parallel processing
num_cores <- detectCores() -1  # Adjust the number of cores as per your system's capabilities

# Register the parallel backend
cl <- makeCluster(num_cores)
registerDoParallel(cl)

# Define the permuteTiming function
pseudo_effects <- foreach(i = 1:1000, 
    .combine = 'cbind', 
    .packages=c( "data.table", "fixest", "ggplot2", 
   "patchwork", "tictoc", "RPushbullet" )) %dopar% {
  res <- permuteTiming()
  return( res )
}
stopCluster(cl)


df1 <- data.frame( t(pseudo_effects) )
names(df1) <- "pseudo_effects"


true_effect = 0.36
f1 = ggplot( df1, aes(x = pseudo_effects)) +
  geom_density() +
  geom_vline(xintercept = true_effect, 
             colour = 'red') +
  xlim( c( NA, 0.4 ) )+
  labs(x = "Placebo effect", y = "Density",
       subtitle = "Placebo Treatment Timing") +
  annotate(geom = "text", x = 0.34, 
           y = 5, label = "Observed effect",
           color = "black", angle = 90, size = 8) + 
  theme( legend.position = 'none', 
         axis.title = element_text( size = 20 ), 
         axis.text = element_text(size = 20), 
         plot.subtitle = element_text(size = 22))
ggsave("appendix_figureA1PanelA.pdf", 
       f1,
       width = 6, height = 5,
       device = cairo_pdf)



# %% Placebo 2: permute scheduled labels
permuteSched = \(){
  # draw scheduled status with empirical probability of scheduled in state
  sch_prob_by_state = above_med[, sch[1], .(cellid, state)][, schprob := mean(V1), state]
  sch_prob_by_state[, pseudo_sch := Map(\(p) rbinom(1, 1, prob = p), schprob), cellid]
  # merge it back and compute treatment effect
  dftemp = merge(above_med, sch_prob_by_state[, .(cellid, pseudo_sch)], by = 'cellid')
  # real treatment timing, fake scheduled status
  dftemp[, Dfake := (year >= first_pesa_exposure) * pseudo_sch]
  m2 = feols(forest_index ~ Dfake | cellid[t] + styear, data = dftemp, cluster = "blk")
  m2$coefficients
}


# Register the parallel backend
cl <- makeCluster(num_cores)
registerDoParallel(cl)
# Define the permuteTiming function
pseudo_effects2 <- foreach(i = 1:1000, .combine = 'cbind', 
                       .packages=c( "data.table", "fixest", "ggplot2", 
                         "patchwork", "tictoc", "RPushbullet" )) %dopar% {
  res <- permuteSched()
  return( res )
}
stopCluster(cl)

# Dataframe of fake results
df2 <- data.frame( t(pseudo_effects2) )
names(df2) <- "pseudo_effects"


true_effect = 0.36
f2 = ggplot( df2, aes(x = pseudo_effects)) +
  geom_density() +
  geom_vline(xintercept = true_effect, colour = 'red') +
  xlim(c(NA, 0.4))+
  labs(x = "Placebo effect", y = "Density",
       subtitle = "Placebo scheduled status") +
  annotate(geom = "text", x = 0.34, 
           y = 5, label = "Observed effect",
           color = "black", angle = 90, size = 8) + 
  theme( legend.position = 'none', 
         axis.title = element_text( size = 20 ), 
         axis.text = element_text(size = 20), 
         plot.subtitle = element_text(size = 22))

ggsave("appendix_figureA1PanelB.pdf", f2,
       width = 6, height = 5,
       device = cairo_pdf)

#------------------------------------------------------------------------------

